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Abstract. 

The critical state problem in type-II superconductivity is described theoretically 
by a direct optimization method, which allows a straightforward treatment for non 
idealized geometries. 

Based on Faraday's law and the principle of minimum entropy production, the 
magnetic history is built up just by a constrained minimization of the field changes 
along the process. Constraints are in the form J £ A, with J the electric current 
density and A some bounded set. This incorporates the vortex pinning and interaction 
phenomena and may be used for the modelling of anisotropy, inhomogeneities and flux 
cutting interactions. 

In this work, our variational statement is posed on the finite element discretization 
and provides a minimal tool for investigating the effects of the sample's topology 
on the field penetration patterns. Simulations of (i) the contraction and splitting of 
boundaries between the flux free and penetrated regions, (ii) the effect of granularity on 
the superconducting properties, (iii) the influence of defects, and (iv) surface curvature 
phenomena are presented. 



1. INTRODUCTION 

Magnetization measurements related to the flux penetration and exit in type-II 
superconductors have been a benchmark against which various models of the interaction 
between flux quanta and the underlying pinning structure are tested(T]. The critical 
state theory [2j provides a quite intuitive and mathematically simple framework for 
analyzing the above mentioned macroscopic experiments. In brief, one is just led to solve 
Ampere's law in the form | V x = J^., with a material dependent parameter, the so- 
called critical current density. This becomes a trivial problem for idealized geometries 
(circular cylinders and slabs in longitudinal applied field), but may be really tough 
when one tries to simulate real world problems. This includes to allow for arbitrary 
sample shapes, spatial inhomogeneities, and more general restrictions on the current 
density than the original critical state equation J = ±Jc, for one-dimensional systems. 
Several numerical methods have been developed, which partially include the previous 
requirements, and may be classified in two groups. 
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On the one side, some models deal with the full electromagnetic problem which 
arises when flux is changed within the superconductor, by means of a definite current- 
voltage law (typically in the form E = Ec{J/Jc)"', n being a phenomenological 
exponent) IH E]- This kind of theory leads to diffusion equations for the 
electromagnetic fields, which have been solved by a number of techniques. They do 
never hold a true critical state, unless for the limit — > oo. 

On the other hand, the general critical state theories avoid using intermediate 
electric fields and just concentrate on the metastable states which can be generated by 
applying external magnetic field changes [HI EI- The material properties are expressed in 
the form of the restriction for the current density vector J G A, with A some bounded 
set, as dictated by the physics of the problem. To be specific, the system either allows 
dissipativeless current flow for J within the region A, or undergoes a transition to a very 
high dissipation regime (infinite in the mathematical limit), which is instantly rejected 
when J ^ A. As it has been discussed elsewhere jH], the selection of A is made on the 
basis of the physical interactions between vortices and the substrate pinning potential. 
Thus, a restriction on J locally perpendicular to the magnetic field (Jc±) relates to the 
pinning interaction of line vortices, and a restriction on the component of J parallel to 
the magnetic field {Jc\\), comes from intervortex tilt interactions. Also, material isotropy 
can be incorporated by choosing A as a circle, or anisotropy, for instance, selecting A 
in the form of an ellipse. Notice that the one-dimensional critical state hypothesis is 
included as the particular case A = [— J^, J^]. The actual solution J = ±7^,0, which 
was brilliantly intuited by C. P. Bean|j2j, arises from the variational interpretation of 
Faraday's law as a rule of minimum magnetic field changes. 

Along the past years, minimum principle formulations for the magnetic behavior 
of superconductors have been exploited, allowing a simple interpretation of many 
experimentals factsfHl El El IHl- These include the analysis of crossed field and rotation 
field experiments, anisotropic responses and the influence of sample's geometry. 

It should be emphasized that generalized variational methods are not merely 
a substitute to differential equation statements of classical physical problems. 
Additionally, they provide a method to generate laws which incorporate singularities, 
discontinuities or inequality constraints, which can be hardly dealt with by local 
differential statements. This is, for instance, the case of type-II superconductors, which 
can be described by the limit n — oo in the power-law current-voltage characteristic 
(multivalued function). Such benefit was already treated in previous issues, where we 
basically discussed a number of physical observations related to the specific set A in the 
theory. As for the sample's geometry, only the slab geometry (one independent variable) 
was treated. In this paper, we exploit another advantage of the variational statement, 
i.e.: the ease to pose the problem for any geometry on the finite element discretization. 
This allows to simulate a variety of experimental observations related to the sample's 
topology. 

The article is organized as follows. In SecEl we pose the physical problem 
of the critical state in two-dimensional geometry and give some general results for 
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our variational statement. A numerical implementation based on the finite element 
discretization is presented in SeclHl This is then applied to several case studies in SecHl 
i.e.: a cylinder with arbitrary cross section, a granular sample, and a sample with a 
point defect. Further applications and possible extensions of the theory are discussed 
in SecEl 



2. STATEMENT OF THE PROBLEM 



In first place, we will define the macroscopic fields involved in the critical state regime 
and their relation to measurable quantities. Our work has been developed within 
the assumption that the sample's response is dominated by the so-called irreversible 
contribution. This means that the equilibrium response of the flux line lattice in the 
absence of pinning centers as well as any possible surface currents are neglected. For 
global magnetization measurements in which the integrated magnetic moment of the 
sample is measured, this is a very good approximation within the region Hd ^ ^ 
Hc2, unless for weak pinning specimens (see Ref.|10j and references therein). Then, the 
coarse-grained electrodynamics if formulated in terms of 

(i) The fiux density B{x) within the sample is the average of the microscopic field 
intensity h, over a volume which encloses a big number of vortices. 

(ii) The magnetic field H{x), in the absence of equilibrium magnetization, is linearly 
connected to B hj B = ^qH . 

— * — * 

(iii) The averaged current density may be calculated from J = V x H. 

Finally, the connection to observable quantities is done by characterizing the external 
field sources and the sample's response. 

(iv) On neglecting finite size effects, the magnetic source [Hs) enters as a boundary 
condition for the flux density at the surface of the sample. The tangential 
component is continuous Bxi^sur f ace) = HqHs^t- 

(v) Then, the measured magnetic moment of the sample per unit volume is M = 
{B)/fio — Hs = (H) — Hs- Here, the average concerns the whole volume of the 
superconductor. 



2.1. The critical state for long samples 

The magnetization of long specimens with arbitrary cross section was already treated by 
Campbell & Evetts^. A number of general properties of the field penetration profiles 
were obtained on the basis of the critical state equation, which according to the upper 
guidelines may be expressed as 

%(f)^^.= (o. 0). (1) 

However, this method is only valid for very simple geometries (f.i.: elliptic and 
rectangular sections) as, in other cases, the free boundary between the critical (| J| = Jc) 
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Figure 1. Superconductor with disconnected subcritical regions established by 
application of an increasing magnetic field. 

and subcritical (J = 0) regions is unknown, and at least, a reasonable guess is needed. 
This may become very tough as the boundary could even split, giving place to multiply 
connected topology. This is depicted in FigQ where two subcritical regions appear for 
the initial magnetization process of a cylinder with concavities in the cross section. In 
general, a free boundary problem will arise between regions either affected or not by 
flux changes in a definite manipulation. Thus, the previous comment extends to the 
fully pentrated regime when the applied excitation is reversed. 

Additionally, the underlying assumption \J\ = Jc or is clearly insufficient for 
describing the full phenomenological scenario. This requires a more general constitutive 
law of the kind J e A|H]. 

Below, we describe the minimal model which accounts for all the aforementioned 
facts. Free boundaries are no longer a shortcoming, as they will arise from the solution 
of our statement, which easily applies for any shape of the sample. Special emphasis 
will be made on this fact through the application to samples with corners, voids, and 
granular properties. 

2.2. Variational statement 

We proceed along the lines established in previous work[Zj. Recall that Faraday's law 
in time discretized form is equivalent to minimizing the spatially integrated magnetic 
field changes over successive time layers, and that the superconducting properties enter 
in the form of the restriction J G A. 

Consider an infinite cylinder of arbitrary cross section subjected to an applied 
magnetic field parallel to the axis. By virtue of the symmetry, the internal field will 
also be parallel to the axis and continuous across the surface, and the induced current 
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density will flow along the basal plane (XY in what follows). Thus, the Critical State 
problem may be formulated as the condition for evolutionary fleld penetration proflles: 

minimize C[i/n+i (■"?)] = J l-f^n+i — -f^n|^ for m G Aj_ . (2) 
In this equation, the following notation is used: 
fl stands for the planar cross section region. 

i^n+i is the inner Z component of the magnetic fleld at the time layer (n + l)6t. 

is the given previous proflle. 

u is the control variable. In this problem u = gradif = {—Jy,Jx)/Jc, i-e.: the 
current density is measured in units of Jc- 

A±_ is the control region (restriction for u, 90° rotation of A). 

In former research, we have solved Eq.((21) in terms of associated Hamilton differential 
equations within the Optimal Control theory. As we restricted our attention to the slab 
geometry in parallel fleld, a system of ordinary differential equations with prescribed 
boundary conditions had to be solved. Here, the situation is a bit more involved, as 
a partial differential equation problem arises. Nevertheless, the same underlying ideas 
may be used. 

2.3. Optimal control formulation 

On using the formalism of Classical Mechanics, we deflne a Hamiltonian density, 
containing the Lagrangian density to be minimized and associated auxiliary momenta 

n{Hn+i,u,p,x) =p-u- ^|i^„+l -if„P . 

Denoting by H*j^i{x), p*{x) and u*{x) the optimal solution functions (i.e. minimizing 
C and satisfying the control system), the optimal control equations are 

OT~C 

grad = — = u , (3) 

div r = = H:_,, - H4x) . (4) 

and the algebraic condition of maximality 

niH,u*,p,x) >niH,u,p,x) VmgAx. 

It is apparent that this condition is satisfled for a control vector u with maximum 
projection over the associated momentum p. This is completely equivalent to the 
problems with one independent variable. Within Q one has u* = p*/p* (or 0), with the 
physical counterpart \J\ = Jc (or 0). In general, one can still identify the critical state 
condition as 

u* G doA^ , 
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with doA± standing for the union of the boundary of the region Aj^ and the zero point 
{doA^ = {dA^}[j{6}). 

In the case of the slab geometry, the main technical difficulty for solving the problem 
was the existence of a free (unknown) point dividing the region penetrated by flux 
changes and the region free of flux variations at a particular step. In the current case 
of interest (more than one independent variable) the free boundary is either a curve 
(2D systems) or even a surface (3D systems). Even more, as shown above (see Fig^, 
the free boundary is quite commonly defined by a collection of disconnected parts. 
Nevertheless, as it was recognized by Prigozhinjn], the use of variational statements 
offers the advantage of applying direct numerical minimization methods in which the 
distinction between both regions is unnecessary. In this paper, we will also make use of 
this benefit for our formulation. We will introduce an optimization method, which 
operates over the whole superconducting region and from which the free boundary 
automatically arises. However, some further comments on the Hamiltonian formulation 
are convenient as they will provide the basis for the physical interpretation of our 
calculations. The boundary conditions to be satisfied by the canonical variables are 
a generalization of the case for one independent variable. On the one side, owing to 
the absence of demagnetizing effects H*_^_^ is determined on the sample surface dQ by 
continuity of the external applied field. 

On the other hand, the free boundary surrounding the regions unaffected by flux 
variations (i/„+i(r) — Hn(T) = 0) is characterized, according to its undetermined 
position by 

r(r) = . 

Notice, that the general situation F = IJ Fj is allowed. 



2.4- Physical interpretation of the momentum vector p 

In this paragraph we show that the momentum p, which has been introduced in the 
optimal control formulation as a Lagrange multiplier, has a clear physical interpretation. 
Recall that, in the absence of flux cutting phenomena, the continuity equation for the 
magnetic fleld may be written as 
dH 

— + diY Jh = , (5) 

where Jh stands for a magnetic fleld current density (not to be mistaken for the electric 
current density!). In simple words, flux variations within a region of the sample are due 
to the entry or exit of flux lines across its boundary. Now, comparing Eq.® and Eq.(jll) 
we conclude 

p = -JnSt , 

with 6t the time step in our discretized model. Thus, solving for p allows to obtain the 
trajectories followed by vortices {flux channels hereafter) in the magnetization process. 
Recall that, isotropy provided, flux channels should be perpendicular to the sample's 
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Figure 2. Flux channels for the vortex penetration in a superconductor when magnetic 
field is increased. 

surface and finish at the free boundary F, which surrounds the regions unaffected by 
flux variations. This point has been illustrated in Fig|21 where we display the integral 
lines for the vector p in the magnetization process of a superconductor with arbitary 
cross section. Dashed style has been used for the so-called d-lines[4\ (qualitative in our 
plot), which are never crossed by the incoming flux tubes. Continuous contours are used 
for the boundary between critical and subcritical regions. More details on the actual 
calculations leading to this plot are given later in the text. 

3. NUMERICAL SOLUTION 

In order to provide a wide application range of this work we describe a triangular 
Finite Element implementation of the numerical procedure for solving Eq.Q. Although 
rectangular meshes are simple and intuitive, triangulation simplifies to deal with any 
geometry for the region Q (sample's cross section). From the practical point of view, a 
big number of software facilities for the triangulation of any specified geometry can be 
found (and downloaded from the web). On the other hand, the fundamental aspects of 
the method may be found in numerous textbooks on the use of Finite Element Methods 
in electromagnetism (see f.i.fHEl)- Here, we will just recall some notational aspects 
and a few relevant properties for our work. 

3.1. Triangular mesh 

As illustrated in FiglHl the 2 dimensional region Q will be meshed by a collection of 
triangles T. A given triangle T G T is defined by three nodes (T = {k, I, m}). The set 
of all nodes is called Af. 
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Figure 3. Triangular mesh defined on a superconducting cylinder with arbitrary cross 
section ft. 



Any physical property given by a function on may be approximated in terms of 
the so-called nodal functions Xk{x,y), which are nothing but the hary centric coordinates 
of the point referred to the triangle T to which this point belongs. On using 

the language of Classical Mechanics, (Afc,Ai,Am,) are the values of normalized masses 
(thus, positive) which one should place at the vertices of T, in order to get the center 
of mass at (x, y) . It is apparent that they provide an affine basis, in terms of which one 
can express the inner points of a given triangle 

Xk,xi,Xm standing for the coordinates of its vertices in some reference frame. The 
approximation of a function on the mesh is now apparent. For instance, the following 
expression may be used for the magnetic field at a given point 

Hix, y) = J2 HkXkix, y) . (6) 

feeT 

Trivially, Hk is nothing but the value of H at the node k. Eq.® will be the starting 
point for obtaining a matrix version of our variational statement (Eq.Q). 

Let us now review some relevant properties of the nodal functions. For a given 
triangle, with area At one has 

Xa = ^^{a^ + Kx + CaV) a = k,l,m 

where 

flfc = Xiym- yi Xm 

h =yi- ym 

Cfc Xm Xl 

and so on, by cyclic permutation. From these expressions one can verify the very 
important property that the nodal functions form a partition of unity, i.e.: 

Afc(f) + Xi{x) + X„Xx) = 1 . 

Based on all this, two fundamental results for our purposes immediately follow 

grad Xa = {ba,Ca)/2AT on T 

and 

^^^-^ "^t/Q if a = (3eT 
Jt "-^^ - \ At/12 if ay^PeT 
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3.2. Optimization problem 

In this paragraph, we show that our variational statement (Eq.Q) may be transformed 
into a nonhnear quadratic optimization problem by using the nodal functions expansion 
of the field (Eq.®) and the above displayed properties of the basis. 

Thus, for the isotropic case, in which the critical current restriction may be written 
l^l = |gradi7| < 1, the magnetization problem consists of finding the set of mesh 
coordinates {Hn+i,k , k G N} for each time layer, so as to 



Minimize 


(-f^n+l,a — H^,a)Iap{Hn+l,p — H^,p) 

T&r a,P£T 


for 


Hn+l,aGapHa+l^l3 < 1 W T E T 

a,PGT 


and 


Hn+l,a = Hn+lfl if « G A/q 


Above, we have used 


lap ~- 




GaP 


= gradAa ■ gradA/3 



and A/q for the subset of boundary nodes, where the magnetic field is determined by the 
boundary conditions. 

Some remarks seem to be convenient before discussing the actual implementation 
of the numerical algorithm stated above. Firstly, we want to emphasize that the method 
is flexible enough so as to incorporate any critical current restriction. Thus, anisotropy 
may be modelled by redefinition of the gradient matrix as 

Gap = (babp + l'^CaCp)/{AAl) , 

with 7 the critical current anisotropy constant (in this case Jcx = Jcy/l)- Spatial 
inhomogeneities in Jc may be dealt with, just by using a triangle dependent scaling 
factor in G^p- Finally, field dependent critical currents would be introduced multiplying 
Gap by some definite function of H at the triangle midpoint (barycenter) . 

From the practical point of view, one will be typically facing a large nonlinear 
and nonlinearly constrained optimization problem. Here, by large it is meant that the 
number of variables, as well as the number of restrictions could be of several thousands. 
Computationally, this seems quite a formidable task for a small machine. Fortunately, 
this kind of problem has been met in the modelling of many physical and scientific 
phenomena (both nature and society tend to optimize and seldom behave as linear 
systems), and one can find very efficient software solutions, even for a personal computer. 

Next, we describe the actual implementation of the previous algorithm which we 
have adopted for the simulations of several physical problems in this work. It is based 
on the computing language MATLAB[13^ and on the fortran package LANCELOT |14j 
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Figure 4. Permanent current density streamlines, obtained for an arbitrarily shaped 
sample by a cycling the external magnetic field. 

for nonlinear optimization problems. We have found that this combination provides a 
rather convenient environment, as one can easily handle the finite element mesh and 
transfer the data structures from one program to the other. The selection of the package 
LANCELOT as opposed to the straightforward use of the MATLAB optimization 
facilities, relates to the higher efficiency of the calculation. Our program is organized 
according to the following steps. 

(I) Definition of the geometry of the region Q. 
(II) Calculation of the finite element triangular mesh on Q. 

(III) Evaluation of the matrix elements Gap and /q/? for the actual mesh. 

(IV) Set up of the standard input file for LANCELOT on the basis of the previous field 
profile Hti{x) and the matrix elements calculated in (III) 

(V) Solution of the optimization problem: get the new profile H^^i{x). 

(VI) Update field profile and iterate from (IV). 

It mus be emphasized that the previous process is performed over the full mesh 
(i.e.: for unknowns i?n+i,k , V A; e and under the constraints \ut\ < 1 V T e T). 
The definition of regions with or without flux variation {Hn+i{x) — H^ix) 7^ or 
H^+i{x) — H^{x) — 0), the boundaries between them (F), and the critical or subcritical 
character = 1 or 0) arise from the optimization process itself. That is to say, the 
totality of mesh points are considered for each iteration, and the optimality process 
determines where to introduce changes or not. This mathematical simulation is just a 
replica of the physical laws ruling out the behavior of the superconductor. Local changes 
are such that global field variations keep minimum, while also fulfilling the minimum 
global entropy production rate (infinite dissipation if J > Jc), and the externally 
imposed boundary conditions. 

4. TOPOLOGICAL EFFECTS 

A rich phenomenology arises in the magnetic field penetration pattern when one 
considers the infiuence of the sample's cross section topology. In this part we report 
on how to understand the concept of critical state for arbitrarily shaped samples. 
First, we consider the influence of the sample's boundary for the so-defined uniform 
superconductors (the current density restriction is homogeneous in space). Secondly, 
we will concentrate on nonuniform superconductors, i.e.: the limitations on the current 
density may change within the sample. 

4-1- Uniform superconductors 
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Figure 5. Field penetration profiles within the superconductor produced by ramping 
up and down of the external field. The vertical axis displays the magnitude of H vs. 
the shape of the sample (horizontal plane) . The external field has been first increased 
from zero, then decreased and increased again. 
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Figure 6. Magnetic hysteresis loop corresponding to the process illustrated in Fig|S| 

4.1.1. General features. The main property which can be assessed by apphcation of 
our optimization technique is that, regardless the shape of the cross section, the flux 
front penetrates through equally spaced contours, which also represent the permanent 
current streamlines. Such a structure is a straightforward consequence of Eq.((H), as the 
model implies a constant gradient, and it has been illustrated in FigEJ Two regions are 
defined by the clockwise or anti-clockwise circulation of the current, as corresponds to 
a magnetic history in which the applied field has been increased and then decreased. 
Notice the boundary line in between. 

Magnetic hysteresis has also been simulated, and one can observe the standard 
Rayleigh-like diamagnetic loop already predicted by CP. Bean for constant Jc in infinite 
slabs or circular cylinders. In fact, the flux penetration and exit that takes place along 
the process (see FiglSl and FiglHI) are nothing but the generalization of the behavior 
for ideal geometries. The external drive pulls up and down across the boundary, while 
maximum gradient profiles are defined within the sample. The interface (free boundary) 
between active and passive regions (as regards flux changes) always corresponds to a 
line of constant distance to the sample's perimeter. 

4.1.2. The effect of corners. A number of outstanding effects related to the penetration 
of magnetic flux in superconductors with corners have been reproduced within the 
framework of our model. We will concentrate on the field profiles for the cross section 
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Figure 7. Splitting of the free boundary in the field penetration process for a sample 
with corners. 

displayed in FigsQ and |21 This example simultaneously holds concave and convex 
corners, and this enables to visualize several pecularities for non ideal geometries. 

FigCI shows a sequence of increasing field penetration profiles. Ten contours are 
plotted, starting with the applied field at the surface. Notice that, as expected, equally 
spaced curves are obtained. However, the particular shape of the outer perimeter leads 
to a splitting process, which takes place as soon as the nearest points of the penetrating 
boundaries contact each other. Further profiles are characterized by a disconnected 
structure in which the free boundary between flux penetrated and flux free regions is 
the union of disjoint curves (F = Fi IJ F2). 

Another group of properties is related to the corners themselves. For the case of 
convex corners, we have observed the formation of the so-called d-lines^ as the locus 
of points where the current paths display sharp bends. Such lines are intuited in Figd 
and have been outlined in FigO Notice that the flux channels do never cross such lines. 
Geometrically, as vortices always move perpendicular to the current streamlines, they 
have to be stopped at the points where J abruptly changes direction. 

For concave corners, a quite different phenomenon is observed. As it was already 
reported for the case of rectangular geometry concave corners become a source of 
magnetic flux in the critical state. This fact is apparent in our simulations. Notice that 
flux channels display divergence at the concavities, because the flux which replenishes a 
finite region must be supplied from a single point. Such behavior is of high importance 
because this point is a natural candidate for the triggering of magnetic instabilities. Just 
recall that the induced electric field at any point may be identified with the amount of 
magnetic flux which has crossed the unit length measured along E in the unit time. 
Thus, E would become infinite for an ideal corner and, in any case, large for a real one. 
As the amount of generated power per unit volume in the critical state is EJc, heat 
removal would first fail at such points (provided the sample is thermally homogeneous) . 

4-2. Nonuniform superconductors 

In this section, we present numerical simulations for samples displaying nonhomogeneous 
behavior of Jc across the section. In principle, any physically meaningful dependence 
Jc{x,y) could be dealt with. Here, we concentrate on two specific cases: granular 
samples and samples with voids. 

4.2.1. Granular samples Granularity has been a main fact in the development of 
superconductivity, mainly since the venue of the high temperature superconducting 
oxides. Owing to their complex chemical structure, these materials have been 
widespread obtained in ceramic form, typically as a collection of superconducting 
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Figure 8. Penetration of magnetic field in a granular sample. The first panel 
illustrates the critical current inhomogeneity (superconducting granules within a poorly 
superconducting matrix). A sequence of field penetration contours follows. 

Figure 9. Penetration of magnetic field in a sample with a void. The magnetic 
field intensity is displayed vertically, grounded on the meshed structure used for the 
simulation. 



grains embedded in a poorly superconducting matrix. Their properties have been 
usually described at a mean field level, characterizing the matrix with some effective 
parameters pn]- However, either analytical or semianalytical statistical theories have 
been developed under certain simplifying assumptions. Here we emphasize that 
the magnetic behavior for any geometrical arrangement may be simulated with our 
numerical approach. The basic idea that magnetic fiux is first screened from the whole 
volume of the sample by means of intergranular currents and further penetrates within 
a collection of basically disconnected grains is visualized. 

FiglHl depicts the successive field penetration profiles for a granular sample with 
rectangular cross section. A collection of cylindrical granules with elliptical cross sections 
has been used, assuming a constant critical current within them (Jcg) which is ten times 
larger than the intergranular limit {Jd = Jcg/10). Notice that, at first, a critical state 
characterized by Jd is built for the intergranular region. The superconducting grains 
are connected in the sense that the field is excluded from the region in between. As soon 
as the intergranular shielding is not intense enough, a double critical state is defined. 
The magnetic fiux penetrates with constant slope Jd for the intergranular region, and 
with the slope Jcg within the grains. 

4-2.2. Samples with voids. We have simulated the field penetration profiles within a 
long superconductor of square cross section with a long cavity of circular shape along the 
axis, and close to one of the edges (see FigE)). The structure of the current distribution 
of such a defect was already discussed by Campbell & EvettsjT] and has been accurately 
reproduced by our optimization principle. Two different simulation methods have been 
used, producing the same results for the magnetic field evolution. 

On the one side, one can describe the whole region as a nonhomogeneous 
superconductor, formed by a circular cylinder with zero critical current Jc,h = {VLh 
will be used for this region), embedded in a matrix with finite critical current Jc^m 7^ 
(f2m will be used for this region). Then, one has Vt = VLh{j Vtm and the set of boundary 
nodes A/o corresponds to the outer square perimeter doVLm- Minimization is performed 
under the restrictions 




V TeT{Qh) 

1 V TeT{nj 
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Notice that equality constraints correspond to the absence of a pinning structure for the 
magnetic flux within the hole. 

On the other hand, one can also solve the problem by excluding the hole from the 
mesh. Then, the boundary of the superconductor splits into an outer boundary do^m 
(square perimeter) and an inner boundary diQm (circular perimeter). Minimization is 
performed according to 

Hn+l,aGaf3Hn+l,l3 < 1 V T 

a, per 

Hn+l,a = Hn+i^fS V tt, /? G C^i^m • 

This means that the magnetic field along the outer boundary is determined by the 
value of the external excitation, and that it must be constant along the boundary of 
the hole. This last condition stems from the fact that the current must flow parallel 
to the surface of the hole. Then gradiJn+i = u is perpendicular to this surface, and 
thus ifn+i(f^i^m) = constant. Notice that this constant is a priori unknown and will be 
determined by the optimization process itself. 

As a further technical detail, before discussing the physics of the problem, we would 
like to comment that the accuracy of the solution may be increased within a reasonable 
computational time, by a specific subdomain decomposition related to the structure of 
the solution itself. Thus, in FigEl we display the actual decomposition of Qm in four 
triangular subregions, according to the well known d-line structure of the square cross 
section, in which sharp bends of the current lines appear along the diagonals. Of course, 
the optimal solution itself should not depend on the specific meshing process, but a much 
larger number of triangles is required for other choices. 

Let us now analyze the physical process of magnetizing such a sample. As it is 
apparent from FigEl the magnetic field levels display a constant gradient within the 
sample. H linearly decreases from a constant value at the surface (the applied field Hs 
which equals -f^n+1,0 for the nth time step), displaying a characteristic plateau at the hole 
area. This plateau changes height along the process, as it is connected to the surface by 
a constant slope region. Thus, the magnetic field within the sample is refresehed in the 
following way. At the beginning, flux penetrates along channels, which are perpendicular 
to the outer perimeter and finish at an equidistant square front. Then, as soon as this 
front reaches the hole a new physical mechanism begins. This has been illustrated 
in Fig^l by means of a vector plot of the field —u = — gradifn+i at a particular 
time step. Recall that, according to the optimal control equations, the magnetic field 
current density Jh is proportional to the vector —u (see EqEland aside). Then, it is 
apparent that, while following the indicated paths, flux lines will preferentially penetrate 
along the highlighted channel. This channel carries all the flux which, on the one side 
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Figure 10. Flux channels in the vicinity of a hole, as represented by a vector plot of 
—u, proportional to the magnetic field current density —Jh (see Eq|Sland aside). A 
collection of vortex lines is depicted within the void. 



replenishes the volume of the hole, and in addition diverges radially through its surface, 
penetrating again into the superconductor. We want to emphasize that, as remarked 
by other authors^l C^I, the magnetic flux outcoming from the hole meets the channels 
coming from the outside boundary in a parabolic d-line. This is a simple consequence of 
the constant critical current ansatz. Discontinuity arises by sharp bending field contours, 
and this takes place at points equidistant to the center of the hole and to a line parallel 
to the horizontal edge of the sample. Trivially, such points are the locus of the parabola 
y = D + R — x^/4i?, with D standing for the vertical coordinate of the center of the 
circle, and R for its radius (the origin has been taken at the center of the square). Notice 
the nice fit of this curve and the place where u changes direction in our simulation. 



5. CONCLUDING REMARKS 



We have applied a direct numerical method for investigating several topological 
properties of the critical state in type-II superconductors. The variational interpretation 
of Faraday's law, constrained by the depinning threshold in the form of the rule J G A 
for the current density, leads to a standard nonlinear optimization problem on the finite 
element discretization. This has been solved by using standard algorithms in the case 
of long cylinders with arbitrary cross section, and the isotropic case (A is a circle). 
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Essentially, the values of the magnetic field at the mesh points (-ffn+i,k) are treated as 
unknowns for each time layer and a minimizer for the global change J^\H^^i — H^\'^ is 
sought. As a result of this method, we obtain the free boundary between the regions 
(un) affected by flux changes, and the circulating current density, which are outputs 
of the theory. The algorithm itself chooses whether i7n+i,k — -f^n,k is zero or not and 
reproduces the critical behavior \J\ = Jc, 0. 

The high flexibility of the theory allows to simulate a number of physical phenomena 
as the influence of arbitrary corners in flux penetration, the splitting of free boundaries, 
spatial inhomogeneities in the critical current, and the calculation of critical state 
properties in multiply connected topologies (samples with voids). 

The fundamental issue \J\ = Jc,0 of our numerical studies for the 2D problem 
confirms the result: gradiJ */Jc = u* = p* /p* , which is a straightforward consequence 
of the optimal control formulation for the variational statement. This reflects a very 
general property {u* G ^oA^), relating to the fact that the cost function does not 
explicitely depend on the control variables. In particular, it can be shown that the still 
challenging full 3D problem is also characterized by the condition \J\ — Jc, 0. This is 
shown below. The 3D critical state problem may be posed as 

Minimize C[Hn+i(x)] = / ^ |#n+i - -^nP 

for V X H^+i = Je in ]R^\Q and V x H^+i e A in Q 

Now, let us assume that the restriction set A is convex (this is not necessary, but 
simplifies the analysis). Then, the optimal control solution V x H*_^i = u* verifies 

c{u*) < c[{i - e)u* + en] v-ueA, o<^<i, 

and this leads to 

lini ^ [C{u* + 9{u - u*)) - C{u*) ] > V e A , 

which one can identify as the condition that the directional derivative for allowed 
variations within A is nonnegative around the minimum. 

If one introduces the actual functional of our problem, the following condition is 

met 

/ ^ {Hn+l{u*) - #n) • (^n+l(^?) - ^n+l(^*)) > ^UEA. 

Then, assuming that a vector field e exists on such that Hn+i{u*) — Hn = —V x e , 
one has 

/ (V X e ) • [Hn+i{u*) - ^„+i(n)] > V u e A . 

On using vector formulas, and assuming that the magnetic field at large distances is 
determined by the external sources, the integral may be transformed to 

/ e-{u*-u)>0 yueA. 
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Eventually, as any candidate control function should verify the external current 
restriction, we get 

/ e-{u*-u)>0 VmgA. 
Jn 

In the case of pointwise controls (the critical current at the point (x,y,z) does not depend 
on sample properties at others) this is equivalent to[T7] 

e-u*>e-u VuGA =^ u* E OqA , 

that is to say, the maximum projection rule implies that the control lies on the boundary 
(in the generalized sense, which admits the trivial solution u* = if allowed by the 
boundary conditions). This ensures that the rule |J| = Jc, is valid for the general 3D 
problem within the isotropic case. 

From the numerical side, the solution of the 3D problem needs a more encompassing 
view point and more sophisticated (though available) mathematical tools. On the 
one side, the finite element mesh must be a tessellation by volumes of some shape. 
A tetrahedral mesh would seem a natural generalization. Then, one can make 
use of tetrahedral nodal functions Aj, which verify analogous properties to their 2D 
counterpart [TTj . In fact, a vectorial built (the so-called Whitney elements) is typically 
used in computational electromagnetism as the right geometrical framework. For 

— * 

instance, the natural basis for the magnetic field {H = J2 HeWe) is obtained from 
We = Afc gradA; — A/gradA^, where e is used for the edge joining nodes k and /. 
Finally, although the current density restriction operates on the finite region Q, one 
must minimize over the whole space IRl This process requires special techniques, such 
as the restriction to some artificial boundary, the use of adaptive meshing, with larger 
and larger elements for distant points, etc. Further development of the theory is planned 
along this line. 
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